A CpG Methylation Signature as a Potential Marker for Early Diagnosis of Hepatocellular Carcinoma From HBV-Related Liver Disease Using Multiplex Bisulfite Sequencing

Background Aberrant methylation of CpG sites served as an epigenetic marker for building diagnostic, prognostic, and recurrence models for hepatocellular carcinoma (HCC). Methods Using Illumina 450K and EPIC Beadchip, we identified 34 CpG sites in peripheral blood mononuclear cell (PBMC) DNA that were differentially methylated in early HCC versus HBV-related liver diseases (HBVLD). We employed multiplex bisulfite sequencing (MBS) based on next-generation sequencing (NGS) to measure methylation of 34 CpG sites in PBMC DNA from 654 patients that were divided into a training set (n = 442) and a test set (n = 212). Using the training set, we selected and built a six-CpG-scorer (namely, cg14171514, cg07721852, cg05166871, cg18087306, cg05213896, and cg18772205), applying least absolute shrinkage and selection operator (LASSO) regression. We performed multivariable analyses of four candidate risk predictors (namely, six-CpG-scorer, age, sex, and AFP level), using 20 times imputation of missing data, non-linearly transformed, and backwards feature selection with logistic regression. The final model’s regression coefficients were calculated according to “Rubin’s Rules”. The diagnostic accuracy of the model was internally validated with a 10,000 bootstrap validation dataset and then applied to the test set for validation. Results The area under the receiver operating characteristic curve (AUROC) of the model was 0.81 (95% CI, 0.77–0.85) and it showed good calibration and decision curve analysis. Using enhanced bootstrap validation, adjusted C-statistics and adjusted Brier score were 0.809 and 0.199, respectively. The model also showed an AUROC value of 0.84 (95% CI 0.79–0.88) of diagnosis for early HCC in the test set. Conclusions Our model based on the six-CpG-scorer was a reliable diagnosis tool for early HCC from HBVLD. The usage of the MBS method can realize large-scale detection of CpG sites in clinical diagnosis of early HCC and benefit the majority of patients.


INTRODUCTION
Hepatocellular carcinoma (HCC) is the most common primary liver tumor with high morbidity and mortality. It has been estimated that hepatitis B virus (HBV) infection causes 50%-80% of HCC cases worldwide (1). The progression of multi-stage hepatocarcinogenesis from chronic HBV infection (CHB) to HBV-related liver cirrhosis (LC) and finally to HBV-related HCC is complex. There is 2%-7% incidence of LC per year in CHB patients (2), and 2%-4% annual incidence of HCC in LC patients (3). Meanwhile, CHB patients can also develop HCC without cirrhosis. The new data show that annual HCC rates ranged 0.03%-1.57% among non-LC Asian CHB male patients (4). In involving 34,952 patients, the study suggested that 2.29% of CHB patients would develop HCC despite hepatitis B surface antigen (HBsAg) seroclearance (5). There are several mechanisms related with HBV-related hepatocarcinogenesis progression including viral regulatory HBV × protein interrupting liver cell proliferation and increasing HBV replication (6), integration of HBV DNA into the host cell genome provoking host cell chromosomal alterations and insertional mutagenesis of cancer genes (7), and host genomic and epigenetic aberrant variation induced by inflammation (8).
The high mortality rate of HCC is mostly due to its discovery at advanced stages. There is an urgent unmet need for biomarkers that can detect early HCC development in CHB and LC liver background. DNA methylation profiles for early stage of HCC in peripheral blood mononuclear cells (PBMC) are different from healthy controls as well as CHB. HCC has a specific DNA methylation signature in easily accessible PBMC and can serve as "noninvasive" biomarkers for detection of early HCC as well as HCC progression (9). Several studies have demonstrated CpG methylation prepared for constructing diagnostic, prognostic, and recurrence models for HCC (10)(11)(12).
However, in these models, measurement methylation of CpGs was performed by pyrosequencing or methylation specific PCR. These methods were laborious and involved fussy work (13) for detection methylation of dozens of CpGs and inevitable increase of inter-batch differences in large samples. This could cause adverse effects on the diagnostic efficacy of the model. Therefore, in this study, we applied multiplex bisulfite sequencing (MBS), which is based on the next-generation sequencing (NGS) method to assess the methylation status of 34 CpGs in 654 samples at the same time. Multivariable methods identified a minimal set of CpGs to achieve optimal prediction. Meanwhile, we handled multiple imputation to complete missing data and bootstrapped training datasets for backward feature selection. The model incorporated both the six-CpG-scorer methylation panel and demographic and clinical characteristics risk factors for early diagnosis onset of HCC from HBV-related liver disease (HBVLD) including CHB and LC.

Data Processing for CpG Selection
The raw IDAT files of Illumina 450K Beadchip data including CHB patients (n = 10) and Barcelona clinic liver cancer staging system (BCLC) 0 (n = 10), BCLC-A (n = 10), BCLC-B (n = 10), and BCLC-C (n = 9) were obtained from Gene Expression Omnibus (accession number: GSE67170) (9), which was prepared to obtain differentially methylated CpGs between CHB and each HCC stage. Forty-four paired PBMC samples were taken from 22 patients before and again after diagnosis of HCC. The 22 patients were diagnosed as LC prior to HCC, and then 5 were diagnosed as HCC BCLC-0 and 17 were BCLC-A. Twenty-two LC and 22 early HCC PBMC samples were subjected to genome-wide DNA methylation assay by using Illumina EPIC Beadchip (data did not published). The R package ChAMP was applied for methylation raw data processing and the paired t tests were used to analyze the differentially methylated CpGs between LC and HCC BCLC-0 and -A. The differentially methylated CpGs with |delta beta|≥ 0.2 and remaining significant after Bonferronicorrected p value < 0.05 were selected.

Patient Study Population
The cross-sectional retrospective study, using the convenient sampling method, was conducted from August 2010 to July 2019 at Beijing You'An Hospital, which is an infectious diseases specialist tertiary hospital in China. A total of 654 patients were included in this study. The demographic and clinical characteristic data of patients with CHB, LC, and early HCC were collected from hospital electronic medical records. The inclusion criterion of CHB was that hepatitis B surface antigen (HBsAg) seropositive status lasted at 6 months or beyond according to the Asia-Pacific clinical practice guidelines on the management of hepatitis B: a 2015 update (14). Diagnosis of hepatitis B-related LC was based on a combination of clinical, laboratory, imaging features, and liver biopsies according to the guideline of prevention and treatment for chronic hepatitis B (2010 Version). Diagnosis of HCC was based on the radiological/ histological criteria according to 2012 EASL clinical practice guidelines: Management of chronic hepatitis B virus infection, Barcelona clinic liver cancer staging system (BCLC) (15). Exclusion criteria were as follows: less than 18 years old; coinfection with human immunodeficiency virus, hepatitis C virus, or hepatitis D virus; coexistence of liver injury caused by drug intake, alcohol consumption, and autoimmune hepatitis; pregnancy; and lactation. The study was approved by the Institutional Ethics Committee of the Beijing You'An Hospital. All participants signed written informed consent on enrolment.

DNA Bisulfite Converted and Multiplex Bisulfite Sequencing
After extraction, PBMC DNA bisulfite was converted with sodium bisulfite according to the protocol of EZ DNA Methylation-Direct Kit (ZYMO Research, Irvine, USA), whose bisulfite conversion rate was >99.5%. Bisulfite-transformed DNA was a single strand and we measured the concentration of the single-strand DNA using a Qubit 2.0 (Thermo, USA, catalog Q10212) to ensure that adequate amounts of single-template DNA were invested for library construction. A panel contained 34 "CpG-free" primer pairs designed to target all candidate CpGs. Library preparation was performed by nest-PCR. Firstround PCR reaction was set up as follows: bisulfite converted DNA 5 ml; forward primer mix (10 mM) 1 ml; reverse primer mix (10 mM) 1 ml; 2×PCR Ready Mix 15 ml (total 25 ml) (KAPA HiFi HotStart ReadyMix PCR). The plate was sealed and PCR was performed in a thermal instrument (BIO-RAD, T100TM) using the following program: 1 cycle of denaturing at 98°C for 3 min, then 27 cycles of denaturing at 98°C for 20 s, annealing at 60°C for 4 min, and final elongation at 72°C for 2 min. Finally, hold at 10°C. The PCR products were checked using electrophoresis in 1% (w/v) agarose gels in TBE buffer (Tris, boric acid, EDTA) stained with ethidium bromide (EB) and visualized under UV light. Then, we used AMPure XP beads to purify the amplicon product. After that, the second-round PCR was performed. PCR reaction was set up as follows: first-round PCR amplicon product (10 ng/ml) 2 ml; universal P7 primer with barcode (10 mM) 1 ml; universal P5 primer with barcode (10 mM) 1 ml; 2×PCR Ready Mix 15 ml (total 30 ml) (Kapa HiFi Ready Mix). The plate was sealed and PCR was performed in a thermal instrument (BIO-RAD, T100TM) using the following program: 1 cycle of denaturing at 98°C for 1 min, then eight cycles of denaturing at 98°C for 20 s, annealing at 60°C for 20 s, elongation at 72°C for 30 s, and a final extension at 72°C for 2 min. Then, we used AMPure XP beads to purify the second-round PCR amplicon product. The libraries were then quantified and pooled. Pairedend sequencing of the library was performed on the HiSeq XTen sequencers (Illumina, San Diego, CA).

Data QC and SNP Calling
MBS raw data were processed by Sangon Biotech. Briefly, raw reads were filtered according to three steps: (a) removing adaptor sequence if reads contains by cutadapt (v 1.

CpG Selection and Assessment
A multi-CpG-based scorer was constructed to diagnose the early HCC patients from HBVLD in the training dataset. The least absolute shrinkage and selection operator (LASSO) with crossvalidation method was applied to select the most significant CpGs from the 34 described in Supplementary Table S1. The process was mainly performed using the "glmnet" package (16) based on the R software (Version 4.0.3) and depicted in Figure 1A.

Statistical Analysis
Firstly, the restricted cubic splines (RCS) were applied to detect the possible nonlinear dependency of the relationship between early HCC and 17 continuous variables, using three knots at the 25th, 50th, and 75th percentiles of the corresponding variable. There were potential threshold associations between early HCC and five indices [age, direct bilirubin (DBil), total protein (TP), albumin (ALB), and hemoglobin] (p-value for nonlinear < 0.001). We have categorized age to the following groups:<43 years and ≥43 years, DBil to <7 mmol/L and ≥7 mmol/L, TP to <65 g/L and ≥65 g/L, ALB to <42 g/L and ≥42 g/L, and hemoglobin to <140 g/L and ≥140 g/L based on the RCS curve. The other four indices [total bilirubin (TBil), g-glutamyltranspeptidase (g-GT), platelet count (PLT), and lymphocyte count (LYM)] were categorized to normal and abnormal groups based on the medical reference value, respectively. The other seven indices [alanine aminotransferase (ALT), aspartate aminotransferase (AST), alkaline phosphatase (ALP), white blood cell count (WBC), monocyte count (MONO), neutrophil count (NUET), and alpha-fetoprotein (AFP)] were logtransformed, respectively (Supplementary Table S1).
Missing values in original dataset were imputed 20 times using imputation with predictive mean matching ("mice" package) (17). For each of the 20 complete dataset, 500 bootstrapping datasets were generated, and backwards feature selection with the Akaike Information Criterion (AIC) was performed on the 10,000 datasets ("rms") (18). Second-order interaction terms were also tested for inclusion. The selected predictors constructed the final model, whose regression coefficients were calculated using "psfmi" package according to "Rubin's Rules" (19,20). The enhanced bootstrap method was to evaluate the stability of the diagnosis model. Discrimination ability was assessed by the area under the receiver operating characteristic curve (AUROC) and C-statistics. Calibration was applied for assessment of agreement between predicted and observed risks across of the population. The nomogram presented the results of the modeling process. Decision curve analysis (DCA) was performed to determine the model clinical usefulness. The construction process referred to this report (21) and the model analysis conforms to the reporting standards of STARD (22) and was depicted in Figure 1B.

Demographic and Characteristics of the Patients
After a strict pathological diagnosis and exclusion process, 168 patients with CHB, 173 patients with liver cirrhosis, 148 patients with HCC BCLC-0 stage, and 165 patients with HCC BCLC-A stage were collected in Beijing You'An hospital and included into this study ( Figure 1A). A total of 654 patients were randomly assigned The analyses were performed on 18 selected predictors measured on 442 parents involved in the training set. Twenty complete datasets were created after 20 times missing imputation, and 17 candidate variables were detected the possible nonlinear dependency of the relationship with early HCC, of which 7 variables were transformed in a non-linear fashion. Resampling 500 times were from each of the 20 complete datasets, leading to a total number of 10,000 bootstrap datasets. A backward feature selection with AIC was repeated on each of the 10,000 bootstrap datasets for selecting the most relevant risk variables for early HCC. The factors chosen at least once during the feature selection procedure constituted the final mode, whose coefficients were estimated using Rubin's Rule from the 20 complete datasets. In internal validation, the model predictiveness and correcting overfit were assessed using 10,000 separate enhanced bootstrapped datasets. The nomogram presented the final mode. Decision curve analysis and clinical impact curves were performed to determine the final model clinical usefulness. The final mode also showed an obvious diagnosis potential in test set (n = 212, 101 early HCC and 111 HBVLD). HBVLD, HBV-related liver disease; HCC, hepatocellular carcinoma; BCLC, Barcelona Clinic Liver Cancer staging system. to the training set (n = 442) and test set (n = 212). Demographic and clinical characteristics of patients are shown in Table 1.

CpG Selection and Panel Signature Building for Early HCC Diagnosis
On the basis of our previous studies (9), we portrayed differentially methylated CpGs between CHB and each of the HCC phases utilizing Limma package with Bonferroni correction. The number of specifically methylated CpGs between CHB and each of the HCC phases included 2,285 for BCLC-0; 2,233 for BCLC-A; 3,345 for BCLC-B; and 23,596 for BCLC-C. There were 326 differentially methylated CpGs that could specifically distinguish early HCC (BCLC-0 and BCLC-A stages) from CHB ( Figure 2A). Moreover, 20 robust CpG  Tables  S2, S3). Thirty-four CpGs were reduced to six potential predictors using the LASSO regression model. The cross-validated error plot and a coefficient profile plot of the LASSO regression model were produced (Figures 2B, C). The logistic regression was used for building the six-CpG-scorer: where risk score = −0.87 -3.73 × cg14171514 + 2.58 × cg07721852 + 6.91 × cg05166871 − 9.85× cg18087306 + 4.50 × cg05213896 + 4.39 × cg18772205, where values of each CpG were methylation ratio measured from MBS. We then applied this formula to calculate the risk score for early HCC of each patient based on their individual six CpG methylation ratio. The risk score was significantly increased in the early HCC samples versus the HBVLD samples (p < 2.22 × 10 −16 ) ( Figure 2D). The six CpGs and their combination six-CpG-scorer also showed diagnostic accuracy ( Figure 2E and Supplementary Table S4). According to determining of maximum Youden index, 0 severed as the optimal cutoff point of risk score. Therefore, we classified those patients with risk  score < 0 as low-risk group, and those with risk score ≥ 0 as highrisk group. The distribution of demographic and clinical characteristics did not vary significantly between the high-risk and low-risk group (Supplementary Table S5).

Six-CpG-Scorer Signature Was Independent of Clinical Factors
The univariate logistic analysis was performed in 10,000 bootstrap datasets. Six of the 18 candidate variables indicated a higher early HCC risk; among these were higher age (≥43 years), male, individuals with higher AFP, higher six-CpG-scorer, lower TP (<65 g/L), and lower TBil ( Table 2). Four variables were selected by the backward feature selection procedure in 10,000 bootstrap datasets ( Table 2). Age, sex, AFP, and six-CpG-scorer signature were independent risk factors for early HCC. The six-CpG-scorer also showed significantly higher predictive accuracy than AFP (Supplementary Table S6) and other demographic and clinical risk factors.

Estimation for the C-Statistics and Brier Score by Bootstrap Validation
Internal validation was performed using the 500 times resampling enhanced bootstrap method from each of the 20 complete datasets. The result showed negligible model optimism. The apparent C-statistics and apparent Brier score was 0.805 and 0.200, respectively. The optimism of the C-statistics and Brier score was −0.0042 and 0.00164, respectively. The adjusted Cstatistics and Brier score was 0.809 and 0.199, respectively.

Diagnostic Performance and Clinical Usefulness of eHCC Nomogram
The AUROC of the model was 0.81 (95% CI, 0.77-0.85) in the training set. The sensitivity, specificity, positive predictive value (PPV), and negative predictive value (NPV) when used in differentiating the early HCC from HBVLD were 70.0%, 77.8%, 74.4%, and 73.6%, respectively. The calibration curve of the eHCC nomogram for the probability of early HCC demonstrated good agreement between prediction and observation in the training set ( Figure 3B). The decision curve analysis of the eHCC nomogram and that for the model without the six-CpG-scorer is presented in Figure 3C. The DCA showed that if the threshold probability of a patient or doctor was >10%, using eHCC nomogram to predict early HCC adds more benefit than either the treat-all-patients scheme or the treat-none scheme ( Figure 3D).

Diagnostic Performance of the eHCC Nomogram in the Test Set
We further enrolled 212 patients including 111 HBVLD and 101 early HCC to serve as the test set for validation of the diagnostic potential. The risk score was significantly increased in the early HCC versus the HBVLD group in the test set (p = 2.7×10 -7 )  Monocyte count variable was truncated with 1st and 99th percentiles and then performed with nonlinear transformation. OR, odds ratio; g-GT, g-glutamyltranspeptidase; ALT, alanine aminotransferase; AST, aspartate aminotransferase. Signif. codes: "***" 0.001, "**" 0.01, "*" 0.05.
( Figure 4A). The eHCC nomogram achieved an AUROC value of 0.84 (95% CI 0.79-0.88) between the early HCC and HBVLD ( Figure 4B). The calibration curve demonstrated good agreement between prediction and observation in early HCC ( Figure 4C). The sensitivity, specificity, PPV, and NPV were used in differentiating the early HCC from HBVLD and were 68.9%, 82.9%, 80.0%, and 71.8%, respectively. The nomogram also indicated good clinical benefits in DCA, which suggested an obvious diagnosis efficacy for early HCC from HBVLD ( Figure 4D).

DISCUSSION
The identification of differentially methylated CpG genome-wide from the Methylation Chip (450K and EPIC) often required a test and/or replication in extra cohorts. The pyrosequencing, bisulfite-conversion-based methylation PCR, PCR cloning, and Sanger sequencing method were laborious and involved fussy work for detection methylation of each CpG site and inevitable increase of inter-batch differences in large samples. The deep sequencing of specific target CpGs often called for MBS, which is based on the NGS method to assess the methylation status at target DNA regions and with very high coverage (13). The important and commendable characteristics of MBS were low cost, low DNA input, and high repeatability and scalability. In addition, the MBS was compatible with common NGS platforms using standard equipment, making it available to most laboratories (23). Construction of amplicon libraries according to a target CpG panel was a key step of MBS. In multiplex PCR using bisulfite converted DNA as template, mixed-base primers, whose pyrimidine base Y contains a "mixture" of cytosines and thymines at the cytosine site and whose purine base R contains a mixture of guanines and adenines at the guanine site, greatly affected the amplification efficiency and led to a decrease in library quality. In the design of this study, we applied 34 "CpGfree" primer pairs (i.e., primers that do not bind to a region containing CpG dinucleotides) targeting each region covering candidate CpGs, which had been considered to produce "nonpreferential" amplification of bisulfite-converted DNA genome (24). We successfully establish a version of the "CpG-free" primer pair panel applied in MBS to simultaneously investigate the methylation status of 34 candidate CpGs across a large sample size (n = 654). The complete and accurate methylation data of 34 candidate CpGs allowed us to screen out and integrate into a six-CpG-scorer with application of LASSO regression. CpG methylation served as an apparatus focused on noninvasive biomarkers for chronic diseases {age-related diabetes (25), diabetic embryopathy (26), coronary heart disease (27), nonalcoholic fatty liver disease (28), or cancer [bladder cancer (29), rectal cancer (30), and HCC (31)]} and had been set up by compelling studies. Our recent study looked at CpG methylation alterations of HBV-related liver disease and developed a model discriminating compensated cirrhosis patients from CHB and decompensated cirrhosis ones based on two CpG biomarkers (32). Our previous study also supported the feasibility of using 350 CpGs for detection of each stages of HCC from healthy control (9). In this study, we successfully established and validated a novel diagnostic nomogram based on the six-CpG-scorer to distinguish early HCC patients from CHB and LC ones. Furthermore, this proposed six-CpG-scorer can diagnose the early HCC better than other demographic and characteristics risk factors. Meanwhile, eHCC nomogram validated considerable diagnostic potential to early HCC and indicated good clinical benefits in DCA.
The CpGs serving as independent risk factors of HCC from mining analysis of Methylation Chip data might point out a novel trace for the exploration of HCC progress mechanism. Based on the Illumina Human Beadchip 27K array, sphingomyelin phosphodiesterase 3 (SMPD3) and neurofilament, heavy polypeptide (NEFH) were found to behave as tumor suppressor genes in HCC after validation in vitro and in vivo (33). Gentilini et al. applied the Methylation 450K BeadChip array and showed that four epigenetically regulated candidate genes [adherens junctions associated protein 1 (AJAP1), adenosine deaminase RNA specific B2 (ADARB2), protein tyrosine phosphatase receptor type N2 (PTPRN2), and sidekick cell adhesion molecule 1 (SDK1)] were potentially involved in the pathogenesis of HCC (34). The finding of AJAP1 was consistent with clinical observation of low AJAP1 expression as an independent factor for predicting disease-free survival (DFS) (35). Mechanically, AJAP1 could block epithelialto-mesenchymal transition (EMT) via suppressing b-catenin/zinc finger E-box binding homeobox 1 (ZEB1) signaling.
In this study, cg14171514 was located within the 5'UTR of the neuroblast differentiation-associated protein (AHNAK) gene, and significantly hypomethylated both in the PBMC DNA of early HCC patients compared with HBVLD ones in this study. Consistently, the aberrant hypomethylation status in the AHNAK gene promoter region from HCC patients' PBMC DNA was verified in methylation-specific PCR (MSP) results of our previous studies (36). Aberrant AHNAK methylation level in PBMC DNA was associated with HCC. In addition, AHNAK mRNA had been reported to be overexpressed in liver cancer tissues by qPCR (36) and the cancer genome atlas (TCGA) data (Supplementary Figure S1). Our novel results from immunoprecipitation in combination with mass spectrometry (IP-MS) analysis strongly indicated that AHNAK protein was involved in and promoted HCC progress (data not published). The cg18087306 was located within the body of Lamin B2 (LMNB2) gene, whose expression level in HCC tissue was higher relative to peritumor tissue. LMNB2 was a promising HCC prognostic and diagnostic biomarker (37).

CONCLUSIONS
We had shown characteristic changes of 34 CpGs in HBVLD and early HCC across a clinical cohort, identified the six-CpG-scorer, and validated their diagnostic efficacy. Thus, we proposed that the six-CpG-scorer-based nomogram is a potential non-invasive diagnostic tool for early HCC from HBVLD, and performed internal validation using the 500 times resampling enhanced bootstrap method from each of the 20 complete datasets to evaluate the stability and then applied it to the test set for validation. The external validation of the monogram will be needed in much larger cohorts from different centers or regions to promote the efficacy and stability to further benefit HCC populations.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Ethics Committee of the Beijing You'An Hospital. The patients/participants provided their written informed consent to participate in this study.